{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# B 12-2 多元线性回归（1）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "研究一个因变量与多个自变量的线性关系。\n",
    "\n",
    "## 案例\n",
    "\n",
    "建立血糖至于胰岛素和生长激素的二元线性回归方程\n",
    "\n",
    "### 数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>blood_sugar</th>\n",
       "      <th>insulin</th>\n",
       "      <th>GH</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>12.21</td>\n",
       "      <td>15.2</td>\n",
       "      <td>9.51</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>14.54</td>\n",
       "      <td>16.7</td>\n",
       "      <td>11.43</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>12.27</td>\n",
       "      <td>11.9</td>\n",
       "      <td>7.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>12.04</td>\n",
       "      <td>14.0</td>\n",
       "      <td>12.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>7.88</td>\n",
       "      <td>19.8</td>\n",
       "      <td>2.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>11.10</td>\n",
       "      <td>16.2</td>\n",
       "      <td>13.52</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>10.43</td>\n",
       "      <td>17.0</td>\n",
       "      <td>10.07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>13.32</td>\n",
       "      <td>10.3</td>\n",
       "      <td>18.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>19.59</td>\n",
       "      <td>5.9</td>\n",
       "      <td>13.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>9.05</td>\n",
       "      <td>18.7</td>\n",
       "      <td>9.63</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>6.44</td>\n",
       "      <td>25.1</td>\n",
       "      <td>5.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>9.49</td>\n",
       "      <td>16.4</td>\n",
       "      <td>4.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>10.16</td>\n",
       "      <td>22.0</td>\n",
       "      <td>2.16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>8.38</td>\n",
       "      <td>23.1</td>\n",
       "      <td>4.26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>8.49</td>\n",
       "      <td>23.2</td>\n",
       "      <td>3.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>7.71</td>\n",
       "      <td>25.0</td>\n",
       "      <td>7.34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>11.38</td>\n",
       "      <td>16.8</td>\n",
       "      <td>12.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>10.82</td>\n",
       "      <td>11.2</td>\n",
       "      <td>10.88</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>12.49</td>\n",
       "      <td>13.7</td>\n",
       "      <td>11.06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>9.21</td>\n",
       "      <td>24.4</td>\n",
       "      <td>9.16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     blood_sugar  insulin     GH\n",
       "IID                             \n",
       "1          12.21     15.2   9.51\n",
       "2          14.54     16.7  11.43\n",
       "3          12.27     11.9   7.53\n",
       "4          12.04     14.0  12.17\n",
       "5           7.88     19.8   2.33\n",
       "6          11.10     16.2  13.52\n",
       "7          10.43     17.0  10.07\n",
       "8          13.32     10.3  18.89\n",
       "9          19.59      5.9  13.14\n",
       "10          9.05     18.7   9.63\n",
       "11          6.44     25.1   5.10\n",
       "12          9.49     16.4   4.53\n",
       "13         10.16     22.0   2.16\n",
       "14          8.38     23.1   4.26\n",
       "15          8.49     23.2   3.42\n",
       "16          7.71     25.0   7.34\n",
       "17         11.38     16.8  12.75\n",
       "18         10.82     11.2  10.88\n",
       "19         12.49     13.7  11.06\n",
       "20          9.21     24.4   9.16"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv(\"B_12_2-data.csv\", index_col=0)\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 线性回归\n",
    "\n",
    "多元线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equation of regression:\n",
      "Y = 17.01 + -0.41 * insulin + 0.10 * GH\n",
      "Fit result:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>blood_sugar</td>   <th>  R-squared:         </th> <td>   0.717</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.684</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   21.54</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Mon, 09 Dec 2024</td> <th>  Prob (F-statistic):</th> <td>2.19e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>17:03:03</td>     <th>  Log-Likelihood:    </th> <td> -36.713</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    20</td>      <th>  AIC:               </th> <td>   79.43</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    17</td>      <th>  BIC:               </th> <td>   82.41</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   17.0108</td> <td>    2.472</td> <td>    6.880</td> <td> 0.000</td> <td>   11.795</td> <td>   22.227</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>insulin</th>   <td>   -0.4059</td> <td>    0.094</td> <td>   -4.313</td> <td> 0.000</td> <td>   -0.604</td> <td>   -0.207</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>GH</th>        <td>    0.0977</td> <td>    0.116</td> <td>    0.843</td> <td> 0.411</td> <td>   -0.147</td> <td>    0.342</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 4.628</td> <th>  Durbin-Watson:     </th> <td>   2.243</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.099</td> <th>  Jarque-Bera (JB):  </th> <td>   2.722</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.868</td> <th>  Prob(JB):          </th> <td>   0.256</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.502</td> <th>  Cond. No.          </th> <td>    134.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:}    &   blood\\_sugar   & \\textbf{  R-squared:         } &     0.717   \\\\\n",
       "\\textbf{Model:}            &       OLS        & \\textbf{  Adj. R-squared:    } &     0.684   \\\\\n",
       "\\textbf{Method:}           &  Least Squares   & \\textbf{  F-statistic:       } &     21.54   \\\\\n",
       "\\textbf{Date:}             & Mon, 09 Dec 2024 & \\textbf{  Prob (F-statistic):} &  2.19e-05   \\\\\n",
       "\\textbf{Time:}             &     17:03:03     & \\textbf{  Log-Likelihood:    } &   -36.713   \\\\\n",
       "\\textbf{No. Observations:} &          20      & \\textbf{  AIC:               } &     79.43   \\\\\n",
       "\\textbf{Df Residuals:}     &          17      & \\textbf{  BIC:               } &     82.41   \\\\\n",
       "\\textbf{Df Model:}         &           2      & \\textbf{                     } &             \\\\\n",
       "\\textbf{Covariance Type:}  &    nonrobust     & \\textbf{                     } &             \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lcccccc}\n",
       "                   & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]}  \\\\\n",
       "\\midrule\n",
       "\\textbf{Intercept} &      17.0108  &        2.472     &     6.880  &         0.000        &       11.795    &       22.227     \\\\\n",
       "\\textbf{insulin}   &      -0.4059  &        0.094     &    -4.313  &         0.000        &       -0.604    &       -0.207     \\\\\n",
       "\\textbf{GH}        &       0.0977  &        0.116     &     0.843  &         0.411        &       -0.147    &        0.342     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lclc}\n",
       "\\textbf{Omnibus:}       &  4.628 & \\textbf{  Durbin-Watson:     } &    2.243  \\\\\n",
       "\\textbf{Prob(Omnibus):} &  0.099 & \\textbf{  Jarque-Bera (JB):  } &    2.722  \\\\\n",
       "\\textbf{Skew:}          &  0.868 & \\textbf{  Prob(JB):          } &    0.256  \\\\\n",
       "\\textbf{Kurtosis:}      &  3.502 & \\textbf{  Cond. No.          } &     134.  \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{OLS Regression Results}\n",
       "\\end{center}\n",
       "\n",
       "Notes: \\newline\n",
       " [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            blood_sugar   R-squared:                       0.717\n",
       "Model:                            OLS   Adj. R-squared:                  0.684\n",
       "Method:                 Least Squares   F-statistic:                     21.54\n",
       "Date:                Mon, 09 Dec 2024   Prob (F-statistic):           2.19e-05\n",
       "Time:                        17:03:03   Log-Likelihood:                -36.713\n",
       "No. Observations:                  20   AIC:                             79.43\n",
       "Df Residuals:                      17   BIC:                             82.41\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     17.0108      2.472      6.880      0.000      11.795      22.227\n",
       "insulin       -0.4059      0.094     -4.313      0.000      -0.604      -0.207\n",
       "GH             0.0977      0.116      0.843      0.411      -0.147       0.342\n",
       "==============================================================================\n",
       "Omnibus:                        4.628   Durbin-Watson:                   2.243\n",
       "Prob(Omnibus):                  0.099   Jarque-Bera (JB):                2.722\n",
       "Skew:                           0.868   Prob(JB):                        0.256\n",
       "Kurtosis:                       3.502   Cond. No.                         134.\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from statsmodels.formula.api import ols\n",
    "\n",
    "formula = 'blood_sugar ~ insulin + GH'\n",
    "\n",
    "model = ols(formula, data=df).fit()\n",
    "\n",
    "params: pd.Series = model.params\n",
    "\n",
    "print(f\"\"\"Equation of regression:\n",
    "Y = {params['Intercept']:.2f} + {params['insulin']:.2f} * insulin + {params['GH']:.2f} * GH\"\"\")\n",
    "print(\"Fit result:\")\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由结果可知，回归直线为\n",
    "\n",
    "$ blood\\_sugar = 17.01 - 0.41 \\times insulin + 0.10 \\times GH $\n",
    "\n",
    "### 假设检验\n",
    "\n",
    "#### 1. 模型检验\n",
    "\n",
    "用方差分析方法检验因变量与自变量之间是否存在线性回归关系。\n",
    "\n",
    "- $ H_0 $ : $ \\beta_1 = \\beta_2 = 0 $\n",
    "- $ H_1 $ : $ \\exists \\ i \\in \\{1,2\\} : \\beta_i \\neq 0 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F-statistic: 21.5389\n",
      "P-value: 0.0000\n",
      "Degrees of freedom: 17.0\n"
     ]
    }
   ],
   "source": [
    "f = model.fvalue\n",
    "f_p = model.f_pvalue\n",
    "dof = model.df_resid\n",
    "\n",
    "print(f\"F-statistic: {f:.4f}\")\n",
    "print(f\"P-value: {f_p:.4f}\")\n",
    "print(f\"Degrees of freedom: {dof}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ p < 0.05 $，接受 $ H_1 $，即所求回归方程有统计学意义。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. 偏回归系数检验\n",
    "\n",
    "使用 F 检验或 t 检验对每个自变量的系数进行假设检验。\n",
    "\n",
    "statsmodels 使用 t 检验。\n",
    "\n",
    "假设\n",
    "\n",
    "- $ H_0 $ : $ \\beta_j = 0,\\ j \\in \\{1,2\\} $\n",
    "- $ H_1 $ : $ \\beta_j \\neq 0,\\ j \\in \\{1,2\\} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==== t values of regression parameters ====\n",
      "Intercept    6.880368\n",
      "insulin     -4.312563\n",
      "GH           0.842833\n",
      "dtype: float64\n",
      "==== p values of regression parameters ====\n",
      "Intercept    0.000003\n",
      "insulin      0.000472\n",
      "GH           0.411024\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(\"==== t values of regression parameters ====\")\n",
    "print(model.tvalues)\n",
    "print(\"==== p values of regression parameters ====\")\n",
    "print(model.pvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ p_{insulin} < 0.05,\\ p_{GH} > 0.05 $ ，可以认为胰岛素对血糖的作用有统计学意义，生长激素对血糖的变化无影响。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 标准化回归系数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    1.665335e-16\n",
       "insulin     -7.434092e-01\n",
       "GH           1.452894e-01\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import zscore\n",
    "\n",
    "df_std = df.select_dtypes(include=[np.number]).apply(zscore)\n",
    "\n",
    "model_std = ols(formula, data=df_std).fit()\n",
    "\n",
    "model_std.params"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
